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Abstract. In the setting of the CE/SE method, a new and unified wall boundary 
treatment for the Navier-Stokes and Euler Equations is proposed. In essence, the shear 
stress exerted on the fluid by a wall is modeled as a source term as a part of local space- 
time flux conservation in the vicinity of a wall boundary. When the fluid is inviscid, the 
source term vanishes and the boundary condition reduces to the usual “slip” condition. 
On the other hand, when the fluid is viscous, the source-term effect is consistent with 
the traditional no-slip condition. Numerical results show' that the new r treatment is 
robust, efficient and accurate for viscous and inviscid flows. 


1 Introduction 

In conventional CFD, the slip condition is applied along the wall for inviscid 
flows, while no-slip condition is applied for viscous flows. In this paper, a unified 
wall boundary treatment for inviscid and viscous flows is proposed in the setting 
of the CE/SE method [1-2]. In essence, the shear stress exerted on the fluid by a 
wall is modeled as a source term as a part of local space-time flux conservation 
in the vicinity of a wall boundary. When the fluid is inviscid, the source term 
vanishes and the boundary condition reduces to the usual “slip” condition. On 
the other hand, when the fluid is viscous, the source-term effect is consistent 
with the traditional no-slip condition. The accuracy, efficiency and robustness 
of the new boundary treatment is demonstrated using numerical examples in 
which reflection of an oblique shock on a wall is calculated for both viscous and 
inviscid flows. 


2 The CE/SE Scheme for 2D Viscous Flows 


The 2D Navier-Stokes Equations in conservation form can be written as 


du in df irn dg irn _ df vm 

dt dx dy dx 


& 9vm p. 


m = 1 , 2 , 3, 4 


The column matrices formed by u m , f irn , 9imy fvm, 9vm , m — 
respectively, are: 
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The inviscid fluxes / im , g tyn are functions of u m , while the viscous fluxes f vm 
and g vm are functions of u m ,chz m /<9x, and du^/dy. Let x\ = x, = y and 
x 3 = t be the coordinates of a Euclidean space E 3 . Then the integral form of 
Eq. (1) is 

(p h m • ds = 0, m = 1, 2, 3, 4 (3) 

JS(V) 

where (i) h m =(/im-/t;m, gim-9vm , Wm), (h) is the boundary of a region V 

in £ 3 , and (iii) ds = d<rn with d <7 and n, respectively, being the area and outward 
unit normal of a surface element on S(V), Note that h m can be decomposed 
into the inviscid and viscous parts, i.e., 

hm ~ him h vrn (4) 


where him — and h vm — (/umi 

In the 2D CE/SE method [2], the space-time mesh used (Fig.l(a)) is con- 
structed from a spatial triangle mesh (Fig. 1 (b)). The unknowns are stored at the 
space-time mesh points (marked by squares and circles in Fig. 1(a)) with their 
spatial projections being the centroids of the triangles depicted in Fig.l(b). Each 
space-time mesh point is associated with three conservation elements (CEs) and 
one solution element (SE). As an example, at point G', the associated CEs are 
cylinders E f F f G f D'EFGD (CE (1 )), A’B’G’F'ABGF ( CE ( 2 >) and C'D'G'B'CD 
GB (CE< 3 >). The SE is the union of four planes A' B f C f D f E f F r , GG"B"B, 
GG f! D"D and GG rf F"F, and their immediate neighborhood. As explained in 
[2], inside each SE (such as that associated with point G'), u m , /i m , gim are 
approximated by first-order Taylor’s expansions in x, y and t with the expansion 
coefficients being the functions of the independent unknowns ( u m )c ( Umx)G' 
and (u my )G f ' By imposing numerical analogues of Eq. (3) over CEs and using 
other conditions, a family of CE/SE Euler solvers, referred to as the Euler a-e- 
a-j. 3 schemes, is constructed in [2]. Since the current Navier-Stokes solver is an 
extension of the Euler a-e-a-fi scheme with additional treatments of the viscous 
terms, only these additional treatments will be discussed in this paper. 

To proceed, note that the third (time) component of each vanishes. Thus, 
in calculating the total viscous flux leaving a CE (such as CE^) through its 
boundary, we only need to evaluate the viscous fluxes leaving the CE through 
its four lateral surfaces ABAlB\ BGB f G\ GFG'F ' and FAF'A'. To expain how 
they are evaluated, as an example, we consider the surface FAF'A r (Fig. 1 (c)). 
Let AS = (S x ,S y ,0) be the unit outward normal vector of FAF f A r multiplied 
by its area. Then 

J h vrn • ds ~ Sxfvm ^( w t7i)q > ( u tixx)q ? 4" 

FAF'A' 
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where Q denotes the centroid of F AF r A ! . Because FAF'A' belongs to the SE 
of point A, we assume that 

(^rn)^ ~ (^m) A H" i^mx) a{FQ %a) "1“ (^my) a{VQ 2/a) H" (^mt) A^tjA (6) 

Note that, as explained in [2], {u mt )A is a function of the independent unknowns 
(urr^A, {umx)A and {u yny ) a- Furthermore, it is assumed that 

~ (^mx)Ai (^mt/)(2 ~ (^) 

Since the above scheme is constructed using a triangular mesh, it is compatible 
with unstructured meshes. In addition, the above scheme has been extended to 
become 3D NS solvers, which however is out of the scope of this paper. 


3 Unified Wall Boundary Treatment 

For simplicity, consider the horizontal wall depicted in Fig. 2(a). No grid point is 
placed on the wall. Instead, a ghost point E, which is the mirror image of point 
G with respect to the wall, is introduced. The flow variables u m and their spatial 
derivatives, u mx and u m2/ (m = 1,2, 3, 4), at point E are obtained from those of 
point G by assuming that, at any time f, the flow r fields below and above DF 
are the mirror images of each other, i.e., we have 

f (^m)^ ~ (^m)<2 5 ^ ~ 1 1 2, 4 f (u mx ) £ = ( 'U'rnx^Q i *77 — 1,2,4 
\ ( ) £ = — (^m)(7 , Tfl = 3 ( (^mx)^ ~ ,777 — 3 


\{ u my)£ — ( u my)(7 i m 2, 4 

| i u my) £ — ( u my)(y , 777 = 3 

Using the above conditions and the concepts of the dual space-time mesh 
discussed in [2], one can calculate the flow variables and their spatial derivatives 
at G for the next time step. The calculation involves an integration over three 
conservation elements CE^,i = 1,2,3 (Fig. 2(b)). However, a special treatment 
for the calculation involving CE^ is needed due to the wall boundary lying 
across CE^. 

Let (i) the viscosity fi be a constant; (ii) the wall be an insulated wall; (iii) 
S be the area of the rectangle DFF'D (iv) Q+ and Q_, respectively, be the 
points immediately above and below the centoid Q of DFF'D'. Then, because 
(i) u = v = 0 at the wall; and (ii) the numerical solution is linear in x, y, and t 
within a SE, the mass, momentum and energy fluxes entering into the fluid in 
the triangular cylinder GDFG'D'F ' from the wall form the row matrix 



1 du 
ReL dy ’ 


4 dv 

3 R e L dy 1 



( 9 ) 
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On the other hand, the same four fluxes entering into the fluid in the triangular 
cylinder EFDE f F f D r from the wall form the row matrix 


9w- — R 



I du 
ReL dy 1 


4 dv \ 

3 R7 L d~y' ) Q 


( 10 ) 


Using the mirror image conditions, we have 

(Sk -(£)*■ ('-scSk-C'-abSk <“> 

By using Eqs. (9-11), it is concluded that the mass, momentum, and energy 
fluxes entering into the fluid in GDEFG f D* E f F r , i.e., CE^, from the wall form 
the row matrix 

9w = 9w+ + 9u. = - 2 S- (o, °’ °) Q+ 


In the current treatment, the surviving x-momentum flux in Eq. (12) is treated as 
a source term in flux balance calculation involving conservation element CE^. 
To calculate the flux g un we need to calculate dujdy . For simple laminar flows 
with enough mesh resolution of the boundary layer, because u = 0 at point Q + , 
and u = (uq + ug')/2 at the midpoint M of GG\ we assume 

{u g + u g ,)I2-0 
\ d y/Q + VM-VQ+ 

Note that, g w = 0 for inviscid flows, and the current treatment becomes the 
usual “slip” condition. As such, the present boundary treatment is suitable for 
inviscid as well as viscous flows. 


4 Numerical Results 

Figure 3 is a schematic of a standard test problem for Navier Stokes solvers. 
When the shock is strong and the incident shock angle is large enough, boundary 
layer separation occurs at the shock impinging point. It is assumed that M**, = 
2.0, R e = 2.96 x 10 5 and 0 = 32.6°. The computational domain is [—1,1.7] x 
[0, 1.2]. Clustered cells near the solid wall are employed to resolve the boundary 
layer. Here 28800 isosceles triangles are used. 

The new solid boundary treatment is applied along the wall. In the left 
and upper boundaries, the flow conditions are fixed according to the incoming 
flows. Non-reflective boundary condition is used in the outlet boundary. The 
numerical results are shown in Fig. 4 for pressure contours, pressure distribution 
and skin friction distribution along the wall. The negative skin friction coefficient 
represents the recirculating flow region. The numerical results, including the 
separation location and length, compare favorably with the experimental data 

[3]- 
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The second example is an inviscid shock reflection problem [4]. This problem 
has exact Euler solution. We apply the same wall boundary condition treatment 
with 1 /R e L — > 0. A uniform mesh of 9600 conservation elements is used. The 
numerical results are shown in Fig. 5. These results agree very well with the 
analytical solutions. 

5 Concluding Remarks 

A unified wall boundary treatment for viscous and inviscid flows has been de- 
scribed. In this treatment, the wall boundary condition for inviscid flows becomes 
a special case of that for viscous flows. Numerical results show that this new 
boundary treatment is robust, efficient, and very accurate for viscous as well 
as inviscid flows. This new treatment could also be easily extended to model 
complex flow physics along wall, e.g., flow permeation through a porous wall. 
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Fig. 1 . A schematic of the CE/SE scheme: (a) a space-time mesh; (b) a triangular 
spatial mesh; (c) calculation of viscous fluxes. 
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(b) 


Fig. 2. A schematic for wall boundary 
treatment: (a) a spatial mesh and (b) CEs 
near the wall boundary. 



Fig. 3. Shock/boundary layer interaction 
problem. 



(a) pressure contours. 



(b) pressure distribution along the walk 


3 r I Skin friction distribution on the plats 

\ Numerical results 

ooo Experimental data 

\ 

2 “ 



(c) skin friction along the wall. 

Fig. 4. Numerical and experimental re- 
sults of the shock/boundary layer interac- 
tion. 



(a) pressure contours. 



(b) pressure coefficients at y = 0.5. 

Fig. 5. Euler solutions of the oblique 
shock reflection problem. 
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